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We show by a numerical procedure that a short-range interaction u induces extended two-particle 
states in a two-dimensional random potential. Our procedure treats the interaction as a perturbation 
and solve Dyson's equation exactly in the subspace of doubly occupied sites. We consider long bars 
of several widths and extract the macroscopic localization and correlation lengths by an scaling 
analysis of the renormalized decay length of the bars. For u = 1, the critical disorder found is 
W c = 9.3 ± 0.2, and the critical exponent v — 2.4 ± 0.5. For two non-interacting particles we do not 
find any transition and the localization length is roughly half the one-particle value, as expected. 



PACS number(s): 71.30, 72.15 Rn, 71.55 Jv 

The interplay of disorder and interactions in elec- 
tronic systems has been studied intensively within the 
last two decades Recent experimental results by 

Kravchenko et al. || have presented strong evidence for a 
metal-insulator transition in two-dimensional (2D) high 
mobility Si metal-oxide-semiconductor field-effect tran- 
sistors, and have generated a great deal of interest in the 
problem of the existence of a metallic state in 2D. These 
results have been confirmed by other workers employing 
different materials and designs Q . All these experiments 
show clear indications that strong electron-electron in- 
teractions partially suppress the quantum interference 
effects responsible for localization. At the same time, 
the scaling theory of localization including the combined 
effects of disorder and interactions predicts that a 2D sys- 
tem may remain metallic even in the limit of zero tem- 
perature Q. 

Direct numerical simulations of the problem are ex- 
tremely difficult and, at present, we have to conform with 
solving the simplest related problem, that of just two in- 
teracting particles (TIP) in a random potential. 

In one-dimensional systems, this TIP problem has at- 
tracted a lot of attention since the original works of 
Dorokhov f| and Shepelyansky ||. The problem has 
been approached from different points of view: by using 
a Thouless type block-scaling picture , by mapping the 
TIP problem onto a random matrix problem [fjj, by di- 
rect numerical approaches based on the time evolution of 
wave packets |||), by transfer matrix methods |fTo|| , by 
Green function approaches |ll],[l2| , by exact diagonaliza- 
tion |l3| and by a decimation method [Q . The previous 
works coincide on the existence of a coherent pair prop- 
agation enhancement. 

In this Letter we perform a numerical calculation of 
the TIP problem in a 2D random potential. Our main 
aim is to establish the existence of extended coherent 
two-particle pairs. Although in principle this is against 
the accepted one-particle scaling picture, it should not be 
very surprising, since 2D is the critical dimension for lo- 
calization. The consideration, for example, of spin-orbit 



coupling effects resulted in a metal- insulator transition 

Our algorithm combines an exact implementation of 
von Oppcn et al. approach fll|| with the scaling pro- 
cedure of MacKinnon and Kramer jl6| for the study of 
critical properties of disordered systems. We consider 
long samples and calculate their decay lengths as a func- 
tion of width and disorder energy. These data support 
the assumptions made in the scaling theory and prove the 
existence of a transition from localized to extended states 
at a finite critical disorder. In adition, to strengthen the 
validity of our calculations, we performed the same anal- 
ysis for two non-interacting particles and obtained the 
expected localization length, i.e., half the value of the 
one-particle localization length. 

Due to the computational effort involved we concen- 
trate in the case of bosons with an on-site interaction, 
although we expect that our main conclusion is equally 
valid for fermions. A small test of the dependence of the 
renormalized decay length with the width of the sam- 
ple clearly shows the same trend for the existence of a 
transition as for bosons. 

We consider a system of length L and width M de- 
scribed by the standard Anderson-Hubbard hamiltonian 
for two spinless particles 

H = t ]T \i,j)(k,j\+t J2 MM 

{i,k},j i,{0,l} 

+ J2 +e j )(i,j\ + U = H + U , (1) 

where i (and j, k, I) labels the L x M sites of a square 
lattice, and e; is the random site energy chosen from a 
box distribution with interval [— W/2, W/2]. {i,k} (and 
{j,l}) indicates that the index k (I) runs over the near- 
est neighbor sites of i (j). The hopping matrix element 
t is taken equal to —1 and the lattice constant equal to 
1, which sets the energy and length scales, respectively. 
We choose an on-site interaction with matrix elements 
(i,j\U\k, I) = u5i_kSj,i8i t j, and use lateral periodic bound- 
ary conditions. 
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To obtain the two-particle decay length £ of the hamil- 
tonian (1) we focus on the two-particle GF 

G^(E-Ho-U)- 1 . (2) 

The full GF satisfies Dyson's equation 

G = Go + G UG , (3) 

where Go is the two-particle GF in the absence of in- 
teractions. The eigenvectors and eigenvalues of the one- 
particle problem are enough to construct Go von 
Oppen et al. [ |j"l| noted that for a local interaction we can 
obtain G very efficiently by projecting onto the subspace 
of doubly occupied sites. This is equivalent to solving 
first the non-interacting case and considering the inter- 
action as a perturbation, only acting on the subspace of 
doubly occupied sites. We will refer with a tilde to the 
matrices restricted to this subspace. Solving Eq. (||) for 
G, and taking into account that U = ul, we obtain 

G=(t-uG r 1 G . (4) 

This expression can be evaluated exactly by inverting 
matrices of range equal to the system size, L x M . These 
matrices are full, all elements are relevant, and their in- 
version cannot be alleviated by matrix transfer methods. 
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FIG. 1. lnTr|G| 2 as a function of I for W = 6, and M=2 
(circles), 4 (squares) and 6 (diamonds). Solid symbols corre- 
spond to u = 1 and empty symbols to u = 0. The straight 
lines are fits from which we calculate the corresponding decay 
lengths. 

Let us call G(m 1; m; m2, 71-2) to the matrix element 
of the GF between an initial (doubly occupied) site of 
coordinates (mi , n\ ) , and a final (doubly occupied) site 
of coordinates (m2,«2)- For a given strip of size L x M 
we calculate the following trace 

lnTr|G(0| 2 Mln]T|G(l,z;U)| 2 ) , (5) 

i,3 



with I < L, and where ( ) denotes an average over the 
disorder realizations. We ensure that L is large enough 
to get a linear exponential decay of the trace as a function 
of I, for any disorder W and width M considered. Once 
we reach the exponential regime, we fit the data in this 
regime to a straight line, whose slope a is related to the 
two-particle decay length through £m = —2/ a. 

In Fig. 1 we show In Tr |G| 2 as a function of the length I 
for M = 2 (circles) , 4 (squares) and 6 (diamonds) . Solid 
symbols correspond to the interacting case and empty 
symbols to u = 0. Each data is obtained by averaging 
over at least 300 configurations. The disorder strength 
is W — 6, which for the interacting case already lies 
in the extended regime, as we will see. We can note 
that a linear exponential decay is well established in all 
cases considered. Strictly speaking our results constitute 
a lower bound for the decay length, but the quality of 
the exponential decay implies that a larger decay length 
would have a very small weight and so would correspond 
to very rare two-particle states. To further check the 
validity of our results, we will apply the same procedure 
to non-interacting particles and we will see that we obtain 
the expected results as compared with transfer matrix 
calculations with very long bars (l^Jl^], which are not 
feasible for the interacting problem. 

Finite-size scaling analysis states that the 

renormalized decay length £m/M is a function of a single 
parameter £/M, 

Cm / M = f(M/0 . (6) 

The scaling parameter £ is the two-particle localization 
length in the localized regime, and the two-particle cor- 
relation length in the extended regime. Equation (^J) im- 
plies that in a log-log plot of (;m/M versus M all data 
should collapse in a common curve when translated by an 
amount ln£(VF) along the horizontal axis. This curve has 
a single branch when there is no transition, while it devel- 
ops two separate branches when a transition is present. 
The main aim of this paper is to discern whether £m /M 
collapses into a single or into a double branched curve. 

In Fig. 2, we show the raw data for £m /M as a function 
of the system width M on a double logarithmic scale for 
different values of the disorder. The on-site energy is 
u = 1 and the disorder energies range between W = 6 and 
15, as indicated in the figure. All data were obtained by 
averaging over a number of disorder realizations ranging 
between 300, for the largest M, and 1000, for M = 2. 
We consider the center of the band (E = 0), and a length 
L = 62. The fact that £,m/M increases with M for small 
values of W, while it decreases for large values of W is a 
clear sign of the presence of a transition. 

A scaling analysis of the data shown in Fig. 2 is de- 
picted in Fig. 3, where we have overlapped all points on 
one curve within the accuracy of the data, by shifting the 
data horizontally by a disorder dependent amount, which 
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is determined by a least-square fit procedure |T(| . Fitting 
the data set for W = 15 and M between 4 and 10 to the 
form £m = £ + A/M we obtain the localization length for 
this disorder £(15) = 2.1 ±0.1, which enables us to estab- 
lish the absolute scale of £(W). As there is practically no 
overlap between the two branches, we have to obtain the 
absolute scale for the upper branch by assuming that the 
scaling parameter diverges symmetrically from above and 
bellow at the transition. The existence of two branches 
is a clear indication of a transition. We will see later 
on how for the non-interacting case we obtain only one 
branch, as expected. 
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FIG. 2. £,m/M versus M on a double logarithmic scale for 
two-interacting particles with u — 1 and for the indicated 
values of the disorder parameter W . 

In order to obtain the critical disorder W c and the crit- 
ical exponent v, we performed an statistical analysis of 
the data in the range 8 < W < 10.5 with the Levenberg- 
Marquardt method for nonlinear least-squares models. 
The most likely fit is determined by minimizing the % 2 
statistic of the fitting function, which we choose to be of 
the form 



= y^A t (W- W c ) l M l/v . 
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The critical disorder found for u = 1 is W c = 9.3 ± 
0.2, and the corresponding critical exponent is equal to 
v = 2.4 ± 0.5. The error bar results mainly from the 
uncertainty in the critical disorder. 

The amount by which we must shift the raw data of 
Fig. 2 to get the universal curve of Fig. 3 gives us the 
scaling parameter £ as a function of disorder for TIP. In 
the inset of Fig. 3 we plot the disorder dependence of 
log £ for u = 1 and the center of the band (E = 0) . 

It should be emphasized that we only demonstrate 
the existence of correlated two-particle extended states. 
Probably, there exist many more uncorrelated localized 



states where the two particles are away from each other. 
Our procedure picks up the longest decay length, at the 
energy considered, which is associated with the most de- 
localized states. 
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FIG. 3. Log-log plot of £m/M as a function of £/M for the 
data represented in Fig. 2. Inset: disorder dependence of the 
scaling parameter £. 

The results for u = are qualitatively different from 
those for u = 1. In the non-interacting case, £,m/M de- 
creases with increasing M for all values of the disorder 
W considered, and so there appears only one branch in 
the scaling procedure. In Fig. 4 we represent £m/M as a 
function of the scaling parameter £ divided by M. The 
upper inset shows the disorder dependence of log£. For 
comparison we also represent the one-particle localiza- 
tion length divided by two £i/2 (solid line), which is 
the expected result for relatively strongly localized sys- 
tems at u = 0. We take for £i the value reported by 
MacKinnon and Kramer ]l6|j . The agreement between 
our results and £i/2 is a positive check of the validity 
of our method of calculation. In the range of validity 
of our results we do not obtain any artificial transition 
for non-interacting particles, as previously reported for 
one particle in 2D due to a different interpretation 
of the raw data. We use MacKinnon and Kramer's |l6| ] 
interpretation, which produces no artificial transition ei- 
ther for the one-particle problem or for our two-particle 
states calculations. As a further check, we have applied 
our method to the well studied ID problem. In the lower 
inset of Fig. 4 we show £ versus disorder on a double 
logarithmic scale for two interacting (solid circles) and 
two non-interacting (empty circles) bosons in ID. The 
straight line corresponds to £i/2, where the one-particle 
localization length £i is taken equal to 105/W 2 . We con- 
sider samples with 500 sites. Our results agree with well 
established previous calculations |p|, ^2|Jhl|l . 

The extension of our results to the case of degenerate 
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electrons, so that can be applied to explain the transi- 
tions found experimentally in Refs. j|] and M and to the 
scaling theory of localization including interactions Jl| , is 
a very difficult problem. In three-dimensional systems, 
fmry [Q argued about the existence of an effective two- 
particle mobility edge that would approach the Fermi en- 
ergy faster than the single-particle mobility edge. Our re- 
sults in 2D could be interpreted as due to the existence of 
a two-particle mobility edge which overcomes the single- 
particle mobility edge at W = 0. Before arguing in favor 
of a transition in 2D degenerate disordered systems, one 
should estimate the lifetime of the coherent pairs, ft is 
not clear how fast our two-particle states would decay to 
lower energy excitations. In the case of degenerate elec- 
trons with long-range Coulomb interactions, Talamantes 
et al. Jn|] reported an increase of the localization length 
with respect to the non-interacting case, but they only 
considered the strongly localized regime. 
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FIG. 4. Log-log plot of £m/M as a function of £/M for 
the non-interacting case. Upper inset: disorder dependence 
of log£ along with the one-particle localization length di- 
vided by two, £i/2 (solid line). Lower inset: £ versus W 
for two interacting (solid circles) and two non-interacting 
(empty circles) bosons in ID. The straight line corresponds 
to fi/2 = 105/2W 2 . 

To summarize, we have calculated numerically the de- 
cay length of TIP in 2D disordered bars of several widths. 
The results are consistent with the assumption of a scal- 
ing hypothesis. Through a scaling analysis of the data 
for u = 1, we proved that there is a localized to ex- 
tended transition at a critical disorder W c « 9.3. The 
critical exponent for the localization length is v w 2.4. 
Our method clearly indicates the existence of delocalized 
states for small disorders, although the values of the crit- 
ical disorder and exponent can appreciably change when 
larger system sizes can be handled. 
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